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P RE FACE 


This  report  represents  work  carried  out  in  the  Spring  of 
1972  to  determine  whether  or  not  it  is  possible  to  track  one 
earth  satellite  from  another,  given  the  spacial  coordinates  of 
the  observing  satellite  and  the  angular  coordinates  of  the 

— rVed  at  sPecified  times  .  The  present  study  is  not  intended 
to  be  exhaustive,  but  rather  is  intended  to  establish  the 
feasibility  of  preliminary  earth  satellite  orbit  determination 
from  a  small  number  of  angle-only  observations. 

The  difficult  and  important  questions  of  detection  and  dis¬ 
crimination  against  the  celestial  background  will  not  be  dis¬ 
cussed.  It  will  be  assumed  that  these  difficulties  can  be 
surmounted  and  that  the  necessary  observations  can  be  made. 

To  make  the  problem  tractable,  we  have  chosen  a  few  simple 
examples  of  orbit  configurations  which  were  studied  in  some 
detail.  The  Planetary  Ephemeris  Program1  (PEP)  has  been 
particularly  useful  in  generating  the  ephemerides  needed  in  the 
method  of  preliminary  orbit  determination.  This  extensive 
computer  program  integrates  the  motion,  and  the  partial 
derivatives  of  the  motion  over  many  orbital  periods  about  a 
non-spherical  Earth;  also,  it  takes  into  account  various  other 
perturbations  and  with  subsequent  least  squares  analyses  leads 


v 


to  a  maximum  likelihood  orbit  determination.  The  object  of 
this  study  is  to  determine  the  osculatina  orbital  elements  from 
limited  amounts  (less  than  one  period)  of  angle-only  data  and  to 
compare  these  with  the  values  assumed  in  the  general  PEP  calcu¬ 
lations  . 

Historically,  the  problem  of  finding  the  orbit  of  a  satellite 
from  a  limited  number  of  angle-only  observations  is  very  interest¬ 
ing.  The  following  brief  account  will  serve  as  an  illustration. 

On  January  1,  1801,  Piazzi  discovered  the  asteroid,  Ceres,  after 
which  he  soon  became  ill  and  had  to  cease  his  observations.  The 
asteroid  thereafter  moved  near  the  sun  where  all  observations 
became  impossible.  It  was  feared  that  the  minor  planet  would  be 
lost.  The  methods  of  Laplace  and  Lagrange  lacked  the  required 
precision  to  locate  such  an  object  months  after  its  first  obser¬ 
vation.  It  remained  for  Gauss  to  work  out  a  method,  later  com- 
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pletely  revised  and  generalized  ,  which  allowed  the  rediscovery 
of  Ceres  on  the  night  of  December  7,  1801  by  DeZach.  It  is 
Gauss'  revised  method  which  we  have  used  to  determine  the  approx¬ 
imate  orbital  elements.  A  detailed  discussion  of  it  is  given  for 
reference  and  for  completeness.  A  computer  proqram  1 ased  on 
this  method  was  written  by  one  of  the  authors  ( DME )  in  order  to 
carry  out  the  necessary  computations. 


It  has  been  assumed  throughout  that  it  is  not  practical  to 
follow  a  given  satellite  over  more  than  a  fraction  of  its  period. 
This  is  a  basic  constraint  of  the  problem;  otherwise  PEP  or  some 
similar  procedure  could  be  applied  directly  without  resort  to 
the  intermediate  step  of  preliminary  orbit  determination.  It  is 
to  be  understood  that  once  a  preliminary  orbit  has  been  deter¬ 
mined,  refined  methods  can  then  be  employed  on  subsequent 
observations  (which  are  not  necessarily  contiguous  in  time  but 
which  can  be  unambiguously  associated  with  that  orbit)  to  arrive 


at  a  precise  orbit. 


INTRODUCTION 


I  . 

The  problem  of  accurately  predictina  the  position  of  a 
body  in  its  orbit  based  on  a  minimum  number  of  observations  is, 
in  general,  insoluble  analytically.  In  many  cases  a  two-body 
approximation  is  unrealistic.  For  example,  to  analyze  the 
motions  within  a  star  cluster,  the  theory  of  General  Pertur¬ 
bations  (numerical  integration)  would  have  to  be  applied.  On 
the  other  hand,  most  Solar  System  problems  like  the  Tarth- 
satellite  problem,  can  be  handled  by  the  techniques  of  Special 
Perturbations.  Here  one  assumes  that  the  force  function  .is  of 
the  form  U  +  P.  where  U  represents  the  two-body  approximation 
and  R  the  disturbing  function.  In  general,  the  contributions 
from  R  must  be  small  compared  to  U.  Unfortunately,  tc  employ 
either  technique  one  needs  either  many  observations  covering 
at  least  one  period  or,  equivalently,  the  elements  of  the  os¬ 
culating  orbit.  Assuming  that  this  is  not  the  case,  one  is 
forced  to  neglect  all  perturbations  and  assume  simple  Keplerian 
motion.  Hopefully  the  two-body  approximation  will  at  least 
yield  enough  precision  to  enable  future  observations  to  be 
made  . 

Six  elements  uniquely  define  an  orbit  and  the  position  of 
the  body  describing  the  orbit.  Three  of  these  define  the 
orientation  in  space,  two  of  them  define  the  size  and  shape  of 
che  orbit  and  the  sixth  will  give  the  position  of  the  body  within 
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the  orbit.  In  the  classical  case  of  a  planet  in  an  elliptical 
orbit  about  the  Sun,  the  elements  are  defined  with  respect  to 
the  ecliptic  and  trie  First  Point  of  Aries.  See  Fiqure  1. 

Briefly,  the  six  elements  are:  the  semi-major  axis,  a;  the 
eccentricity,  e;  the  inclination  of  the  plane  of  the  orbit  to 
the  plane  of  the  ecliptic,  i;  the  angle  that  the  major  axis 
makes  with  the  line  of  nodes  (the  longitude  of  perihelion),  ; 
the  angle  that  the  line  of  nodes  makes  with  the  line  from  the 
Sun  to  the  vernal  equinox,  (the  longitude  of  the  ascending 
node),  fi;  and  the  time  of  perihelion  passage,  t.  Six  elements 
have  to  be  found;  hence,  three  observations  of  the  body's  right 
ascension,  and  declination  at  three  different  times  constitute 
the  minimum  data  for  orbit  determination. 

From  this  preliminary  orbit  it  is  possible  to  compile  a 
table  of  predicted  positions,  an  ephemeris,  to  be  used  for 
tracking  the  body  so  that  future  observations  may  be  made. 
Additional  observations  would  then  be  used  presumedly,  to 
improve  the  crbit.  There  are  many  ways  to  accomplish  this  and 
only  the  basic  ideas  will  be  presented  here. 

Assume  the  elements  of  the  preliminary  orbit,  c^,  i  =  1,...6 
have  been  calculated  for  a  geocentric  satellite.  Because  all 
possible  perturbations  were  neglected,  these  elements  are  not 
the  elements  of  the  actual  orbit.  Hence,  any  observed  quantity 
i  where 
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ECLIPTIC 


ORBITAL 

PLANE 


Lij?*  L  Tlie  orientation  of  an  oroit  in  space.  In  the  earth  satellite  case 
the  plane  of  the  earth's  equator  replaces  the  ecliptic. 


=  >*>  (oi,qi,t)  i  =  1,  ...  6 


and  the  qVs  are  the  six  elements  of  the  Earth's  orbit,  will 
differ  from  the  predicted  quantities  ;  at  a  given  time.  The 
;'s  may  be  observations  of  the  nqht  ascension  and  declination 
( -i ,  6 )  ,  the  geocentric  co-ordinates  (x,v,z)  or  the  velocity 
components  (x,y,z),  etc.  Assuming  the  preliminary  orbit  was 
reasonably  close  to  the  actual  one,  the  change  in  the  -  '  xtal 
elements  Aa^  will  be  slight.  The  change  in  .  will  be 


At  =  l 


.  o  O  ■ 

1=1  1 


a  c  .  . 


Then  if  we  let  the  difference  between  the  observed  quantity  and 
the  predicted  quantity  be 


At  =  I  -  t 

o  vt 


we  have. 


where  the  subscript  j  means  the  quantity  was  observed  or  pre¬ 
dicted  at  time  t .  . 
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If  n  =  6,  the  problem  is  a  straightforward  one  of  six 
equations  in  six  unknowns.  If  n  >  6,  then  the  equations  can  be 
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solved  by  least  squares  or  a  similar  technique  In  either  case 
the  Ack  can  then  be  added  to  the  to  yield  "improved"  values 
of  the  orbital  elements.  Theoretically,  the  process  can  be 
carried  out  for  all  sets  of  future  observations.  If  and  only  i 
the  osculating  elements  are  reasonably  accurate  can  we  hope  to 
continually  improve  the  orbital  predictions. 


II. 


THE  METHOD  OF  GAUSS 


Several  methods  have  been  developed  to  solve  the  problem 
of  preliminary  orbit  determination.  Needless  to  say,  until 
the  advent  of  high  speed  computers  the  orbits  so  determined 
were  hardly  to  be  considered  preliminary.  The  basic  methods 
of  Laplace,  of  Lagrange  and  of  Gauss  are  probably  the  most 
useful.  The  first  two  methods  have  been  neglected  in  favor  of 
Gauss'  technique  since  the  latter  makes  no  assumptions  about 
the  time  between  observations  and  readily  lends  itself  to  an 
iterative  procedure. 

Let  (x,y,z)  denote  a  geocentric  rectangular  equatorial 
coordinate  system  as  in  Figure  2:  the  x-axis  is  directed  towards 
the  First  Point  of  Aries.  Let  (f,n,r)  be  a  parallel  coordinate 
system  centered  on  the  Hunter  satellite,  H,  then  the  unit 
vector  u  defining  the  Hunted  satellite  in  the  Hunter  satellite 
reference  system  is 

u  =  (cosacos5)  i  +  (sinacos6)  j  +  sinfk 

where  i,  j,  and  k  are  unit  vectors  along  the  , n ,  axes  re¬ 
spectively  and  ( , 6  ^ )  ,  i  =  1,  3  are  the  three  observations  of 
the  observed  satellite,  h,  from  the  observing  satellite,  H. 

From  Figure  2  we  see  that 

ri  “  »iui  -  R  n> 
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- -  .. 


Fig.  2.  The  vector  relationships  resulting  from  a  hunter  satellite 
centered  coordinate  syste1  1. 
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where  R  =  (X,Y,Z)  is  the  position  of  the  Earth  as  seen  from  the 
Hunter  satellite  and  p  is  the  distance  to  the  observed  satellite. 
We  muse  assume  that  the  equatorial  coordinates  of  the  Earth  are 
known  accurately  at  each  of  the  three  observation  times. 

Since  we  are  assuming  Keplerian  motion  we  know  that  the 
motion  takes  place  in  a  plane;  hence,  one  of  the  unknown  radius 
vectors  may  be  written  as  a  linear  combination  of  the  other  two. 
We  have 

r  2  =  Clrl  +  C3r3  (2) 


where  the  coefficients  c^  c^  are  the  so  called  triangle  ratios. 
If  c^  and  c^  are  known,  then  substitution  of  equation  (1)  into 
equation  (2)  furnishes,  in  component  form,  three  linearly 
independent  equations  in  the  unknowns  p^,  P2'  and  p3'  the  1'-unter_ 
hunted  satellite  distances  at  the  times  of  the  observations. 

Since  the  coefficients  are  unknown,  we  approximate  then  iterate. 
First  we  note  that  the  magnitude  of  the  vector  cross  product  of 
any  two  of  these  radius  vectors  is  double  the  area  of  the  tri¬ 
angle  formed  by  the  Earth  and  the  two  positions  of  the  Hunted 
satellite.  See  Figure  3.  Hence,  we  can  write  the  coefficients 

A 

as  ratios  of  triangle  areas.  Taking  crossed  with  equation 
(2) ,  we  have 


rl  *  r 2  C3rl*  r 3 


c3  =  [r1, r2J/[r1,r3] 


(3) 
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where  the  brackets  denote  triangle  areas.  Similarly  by  cal¬ 
culating  r2*r3  we  find  for  the  coefficient  c.  : 

C1  =  [r2'r3]/[ri'r3] •  <4) 


Let  (r^r^)  denote  the  area  of  the  sector  of  the  ellipse 
bounded  by  the  radius  vectors  r ^ , r ^  and  denote  the  sector- 
triangle  ratios  by 


-  _  (r2'r3} 

Yl  ~  tr2,r3] 


(rl'r3} 
y2  [r1(rJ 


(rl,r2) 

7 3  “  [r .  , r  ]  ' 


(5) 


Then  following  Gauss  we  can  rewrite  equations  (3)  and  (4)  as 
sector-triangle  ratios: 


(r2,r3)y2 

(rl'r3)yl 


(rl,r2)y2 

<rl'r3)y3 


or,  recalling  Kepler's  second  law,  namely  that  the  area  described 
by  the  radius  vector  is  proportional  to  the  time, 


°i  = 


(t3  fc2)y2 
(t3  -  Vyl 


(t2  ~  tl)y2 
(t3  -  Vy3 


(6) 


Encke  developed  a  series  solution  for  the  sector- triangle 

ratios;  his  expansion  to  the  third  order  is, 

-  =  i  x  4m  88m2  8m 2  ,  5312m3  ,  512m22  ,  64m22 

y  3  45  5  945  105  35  *  *  * 

where , 
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m  =  nsec  y 


l  = 


(1  -  COSy) 
2  cos  y 


n  =  uk 


2(tj  ~  ti) 
(ri  +  rj} 


_  (ri‘rj} 

cos  (f4  f^)  -  (r.r.) 


3 


i  3 


secy  = 


r  .  +  r  . 
i  _1 


2  (rir  j  )  1^2cos  (  (f  ^  -  f^/2)) 


k  is  Gauss"  constant,  and  f  is  the  true  anomaly  (the  angle 
swept  out  by  r  since  the  time  of  perigee) .  If  we  now  let 

T ]_  =  k  (t^  ~  1 2 ) 

T2  =  k  ~  fcl^ 
t3  =  k  (t2  -  tx) 


denote  the  time  intervals  between  observations  (in  some  appro¬ 
priate  units  determined  by  k)  ,  we  can  express  the  coefficients 

in  (6)  as, 


to  the  first  order  in  m.  The  coefficients  a^b^a-j,  and  b3 
can  be  calculated  immediately  from  the  input  data.  Substituting 


See  Appendix  A 
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equation  (1)  into  equation  (2)  we  have 

■  2u2  —  C1  2U1  ~  R]_  ^  R2  C3  3U3  —  R3  ^  ( 8 ) 

then,  substituting  the  expressions  in  (7)  for  the  coefficients 
cl'  c3 ' 


Operating  on  this  equation  with  (.u^u^  yields,  after  a  bit  of 
arithmetic,  an  equation  of  the  form: 


p  2  =  A  + 


where  , 

al  [Ri'uix  u3]-[R2-u1xu3]+a3[R3.  Ulxu3] 
[u1-u2xu3J 

bl  [Ri 'uixu3 J +b3 [r3*  ui^  u3 ] 

B  =  - - - - . 

[ul*  u2  U3] 


(9) 


(10) 


A  second  equation  relating  p2  and  r2  is  evident  from  Figure  2. 
By  the  cosine  law 


r 


2 

2 


r: 


2 ( u 2  * R2) . 


(11) 
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We  now  have  two  equations,  (9)  and  (11),  in  two  unknowns  which 

can,  theoretically,  be  solved  for  p2  and  r 2  .  Any  iterative 

procedure  can  be  used  in  the  solution  provided  a  reasonable 

initial  guess  at  r 2  can  be  made.  Alternatively,  one  could 

form  the  following  eighth  order  equation  and  then  use  Eairstow's 
4 

Method  or  an  equivalent  technique  to  find  the  roots.  Sibstituting 
equation  (9)  for  p2  into  equation  (11)  yields, 

r2  (C  +  AD  +  A2)  *2  ~  (BD  +  2AB)  *2  ~  B2  =  0 

where  C  =  R2  and  D  =  -2(u2*R2).  If  the  cot  icients  c1  and  c3 

are  exact,  then  this  equation  will  have  a  trivial  root  r  =  R 

2  2 

at  t  2  -  0.  In  any  event,  one  of  the  roots  will  approximate  the 
trivial  solution  and  either  four  or  six  of  the  roots  will  be 
complex.  The  latter  case  yields  no  solution,  since  the  remain¬ 
ing  real  root  will  be  negative.  In  the  first  case,  then,  five 
of  the  roots  can  be  immediately  disregarded.  It  can  also  be 
shown  by  considering  the  number  of  variations  in  sign  of  the 
coefficients  that  of  the  remaining  three  real  roots  we  must  have 
either  two  positive  and  one  negative  or  all  three  negative. 

Again,  the  case  of  three  negative  roots  will  yield  no  solution. 
Hence,  there  are  two  'possible'  roots  to  this  equation  to  be 
considered.  In  a  well  defined  case  (i.e.,  when  Keplerian  motion 
is  a  good  first  approximation)  one  of  these  possible  roots  will 
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be  outside  of  any  reasonable  limits  to  the  problem.  In  the  Earth 
satellite  case,  for  example,  obvious  limits  would  be  the  radius 
of  the  Earth  and  1  A.U.  (astronomical  unit).  Hence,  in  a  well 
defined  case  there  is  but  one  solution  and  one  orbit  for  the 
satellite  being  observed. 

When  a  solution  has  been  found,  the  coefficients 
c  and  c3  can  be  evaluated  from  equations  (7).  Then  operating 
on  equation  (8)  first  with  (*u2  u3^  ant^  then  with  we 

find  for  the  other  two  Hunter-Hunted  distances: 

C]  [R1  -u^xu-j]  -  [R2-u2xu31  +  c3[R3-u2xu31 
f  1  —  ~  A  ” 

cl[ui‘u2<u3] 

c-j  [R-| -u-jXu^  -  [Rj'U^Uj]  +  c3  [R3  ’ ulx  u2 1 

P  *5  *"  ~  ~  ^ 

c 3  u2; 

Equation  (1)  will  then  give  the  two  radius  vectors  *1  and  r3> 

Because  the  coefficients  c-^  and  c3  were  determined  from 
truncated  series  expressions,  the  solution  p2,  is  only 
approximate.  Hence,  the  solution  might  be  improved  by  iterating 
on  the  coefficients.  First  the  observation  times  and  the  time 
intervals  should  be  corrected  for  light  time  by  using  the  relation 
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J 

1 


t[  =  ti  ~  -577  x  10-^  (days) 


where  o i  is  the  distance  (A.U.)  calculated  in  the  previous 
iteration.  The  sector-triangle  ratios  are  calculated  again 
using  Encke's  series  solution  including  second  and  perhaps, 
third  order  terms  and  using  the  previously  calculated  radius 
vectors  r^  r2  and  r^.  The  "improved"  coefficients  will  be: 


C1  = 


Tly2 

T2yl 


C3 


T  3y2 
T2y3 


If  these  coefficients  agree  to  within  the  desired  accuracy  with 
the  previously  calculated  values,  then  the  iteration  may  be 
stopped.  If  not,  then  the  coefficients  in  equation  (7)  should 
be  recalculated  according  to: 


al  =  T i/x  2 


a3  T  3^T  2 


bi  =  S  \r. 


1  r; 
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Equation  (9)  is  reformed  and  a  second  solution  with  equation  (11) 
is  carried  out.  If  the  three  observations  are  spread  over  a 
relatively  small  section  of  arc,  then  three  iterations  will 
probably  prove  satisfactory.  With  the  final  values  of  and 
the  radius  vectors  r ^  and  r^  are  calculated  by  ecruation  (1)  and 
with  these  the  orbital  elements  are  calculated. 
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Ill . 


CALCULATION  OF  THE  ELEMENTS 


Only  two  radius  vectors  and  the  times  corresponding  to 
these  positions  are  needed  to  calculate  the  six  orbital  elements. 
In  the  heliocentric  problem  it  is  generally  the  ecliptic  elements 
that  are  required.  In  that  case,  the  radius  vectors  r^  = 
(x.,y.,z.)  can  be  transformed  by: 

x "  =  x 

y'  =  ycose  +  zsine 
z'  =  -ysinr  +  zcose 

where  e  is  the  mean  obliquity  of  the  ecliptic.  In  the  geocentric 
problem  the  equator  is  considered  the  fundamental  reference 
plane;  hence  no  transformation  is  required. 

Let  r^  =  (x^y^z-^)  and  r^  =  (x^,y^,z3);  then  the  equation 
of  the  plane  of  the  orbit  will  be 

ax  +  by  +  z  =  0 


where , 


=  ylz3  -  y3zl  _  zlx3  ~  Z3X1 

3  xly3  ‘  x3y1  xly3  '  x3yl 

and  the  equation  of  the  line  of  nodes ,  the  intersection  of  the 
orbital  plane  and  the  plane  of  the  equator  will  be 
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ax  +  by  =  0 . 


But  the  slope  of  the  line  of  nodes  is  just 


tan 


a  _  ylZ3  "  y3Zl 
b  xxz3  -  x3z1 


where  is  the  longitude  of  the  ascenainq  node.  See  Figure  1. 
The  inclination  of  the  orbital  plane  to  the  plane  of  the  equator 
can  be  gotten  from 


cosi 


[  1  + 


a2  + 


b2J1/2 


To  determine  the  quadrants  of  u  and  i,  first  determine  if  the 
points  (x1,y1)  and  (x3>y3)  are  in  different  quadrants.  If  so, 
then  the  direction  of  motion  (direct  or  retrograde)  is  apparent. 
If  the  points  are  in  the  same  quadrant,  then  when  y3/x3  >  Y]/Xq 
the  motion  is  direct  (i  <  90°);  conversely  the  motion  is 
indirect  (i  90°).  The  slope  of  the  line  of  nodes  will  in¬ 

dicate  throuqh  which  quadrants  the  line  passes;  a  positive  slope 
indicates  quadrants  1  and  3,  while  a  negative  slope  indicates 
quadrants  2  and  4.  To  determine  which  of  the  two  quadrants 
contains  the  ascending  node,  it  is  necessary  to  consider  the 
quadrant  and  slope  of  the  third  radius  vector  with  respect  to 
the  line  of  nodes.  If  the  motion  is  direct  and  Y3/x3  tan  > 
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the  satellite  will  pass  through  the  ascending  node  next.  If  the 
motion  is  indirect  and  y3/x3  •’■s  9reater  than  tanft,  then  the 
satellite  will  pass  through  the  descending  node  before  the 
ascending  node.  Thus,  and  i  are  determined.  These  two 
elements  orient  the  orbit  in  space. 

Using  Kepler's  law  of  areas,  we  have  that  the  sector 
area  is, 

2(rJL,r3)  =  kp1/2(t3  -  tx) 

where  p  is  the  semi-lattice  rectum,  p  =  a(l  -  e  )  and 

1/2 

kp  / 2  is  the  magnitude  of  the  area  velocity.  See  Figure  4. 

For  the  triangle  area  we  have 


2  t  ri  ^ r  3  ]  =  r !  r-,sin  (f  0  -  fn) 


1  3‘ 


where  f  is  the  true  anomaly.  Then  using  the  sector  triangle 
ratio  y ^  previously  calculated,  we  can  solve  for  p: 


(r ,  ,  r3) 
2  ~  [rl'r3] 


= 


p 


i/2  _  rlr3Sln(f3  fl)v2 
k(t3  -  t^) 


To  find  the  eccentricity,  e,  we  use  the  equation  of  a 
conic  section,  namely, 
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1  +  ecosf 


Setting  =  p/r^  -  1,  we  form  the  first  of  two  equations  in 
two  unknowns,  e  and  f ^ : 


Similarly , 


ecosf^  =  q^.  (12) 


ecosf3  =  q^  =  p/r^  -  1. 


Since  we  can  calculate  the  difference  in  anomalies  from 
cos(f3  -  f^)  =  (r^*r^)/(r1r^) ,  we  could  rewrite  f^  as 
f^+(f^  -  f^)  and  calculate  esinf^  as  follows 


q^  =  ecos f ^cos ( f ^  -  f^) -esinf^sin  (f^  -  f^) 
q  cos (f  -  f  )-q 

esinf,  =  —7 — 7-p ?  (13) 

1  sin(f3  -  f1) 

Equations  (12)  and  (13)  can  then  be  solved  simultaneously  to 
yield  e  and  f ^ . 

Havinq  found  the  semi-lattice  rectum  and  the  eccentricity 
we  readily  calculate  the  semi-major  axis,  a  from  the  elliptical 


relation : 


The  mean  motion,  n  follows  directly  from. 


1/2 

where  ku  is  just  G,  the  constant  of  gravitation  times  the  sum 
of  the  mass  of  the  Earth  and  the  satellite. 

Next  we  can  find  the  argument  of  latitude  u^  =  u.  +  f 
from  the  equations 

sinu^  =  z^/r^sini 

cosu-^  =  ( x  1  cos  +  y1sinfi)/r1> 

Then  the  argument  of  perigee  is  readily  obtained  from 
to  =  (u^  -  f  )  . 

One  element,  the  time  of  periqee  passage  remains  to  be 
found.  If  the  mean  anomaly,  the  angle  swept  out  by  a  radius 
vector  with  mean  angular  velocity  n  in  the  time  interval  (t  -  t) 
were  known,  we  could  calculate  t  from  the  equation, 

M1  =  ”  (fci  “  1  )  • 

In  order  to  calculate  the  mean  anomaly,  the  eccentric  anomaly 
must  first  be  found.  The  true  anomaly,  already  calculated  for 
time  t.  and  the  eccentric  anomaly, E  are  related  by. 


tan(E1/2)  =  tan(f1/2)|(l  -  e)/(l  +  e)J1,/2. 

In  addition,  from  Figure  4  we  see  that, 

r^cosf^  =  acosE1  -  ae. 

The  second  equation  can  be  used  to  determine  the  quadrant  of 
El*  Finally,  Kepler's  equation  relates  the  mean  and  the 
eccentric  anomalies, 

M1  =  Ei  -  esinE1 . 
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IV.  DISCUSSION  OF  RESULTS  &  CONCLUSIONS 


Two  basic  positions  for  the  hunter  satellite  were  assumed, 
one  times  (42,  164  km)  and  1/4  times  (10,  541  km)  synchronous 
satellite  altitude.  The  following  orbital  configurations  were 
investigated : 

Hunter  Satellite  Hunted  Satellite 


(a) 

1 

sync 

1/4 

sync 

(b) 

1/4 

sync 

1 

sync 

(c) 

1 

sync 

3 

sync 

The  two  satellites  were  placed  in  various  relative 
positions  to  each  other  in  their  respective  orbits.  PEP  was 
then  employed  to  compute  the  earth  centered  inertial  coordinates 
of  each  satellite  and  the  right  ascension  (R.  A.)  and  declination 
(Decl.)  look  angles  from  the  hunter  to  the  hunted  satellite  all 
as  a  function  of  time.  The  quantities  needed  as  input  to  the 
Gauss  program  are  as  follows: 

Ephemeris  Time  R. A.  Dec 1 .  X  Y  Z 


The  X,  Y,  Z  are  the  hunter  centered  inertial  coordinates 
of  the  Earth,  just  the  negatives  of  the  values  computed  by  PEP, 
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which  are  assumed  known  for  all  times.  In  general,  interpolation 
of  the  tabular  values  was  necessary. 

In  the  classic  problem  of  Ceres,  the  planet  moved  a  total 
of  about  3°  in  1  1/3  months  of  observation.  In  configuration 
(c)  the  time  between  observations  was  varied  from  1  hr.  (3°)  to 
6  hr.  (18°).  In  general,  the  method  failed  in  cases  where  the 
total  arc  length  between  r^  and  r^  was  greater  than  about  10°. 
Furthermore,  in  many  configuration  (a)  situations  there  resulted 
two  nearly  equal  positive  roots  of  the  equation  in  r^ /  of  which 
one  leads  to  approximately  correct  orbital  elements,  but  the 
other  one  does  not.  It  was  difficult  to  provide  unambiguous 
instructions  which  told  the  computer  which  one  of  the  two 
roots  was  the  "correct"  one  to  use.  Many  configuration  (a) 
situations  either  failed  outright  or  were  very  slow  to  converge. 

At  first  the  poor  success  rate  in  configuration  (a)  was  thought 
to  be  caused  by  lack  of  coordinate  precision  resulting  in  the 
use  of  linear  interpolation.  However,  when  a  more  accurate 
Hermite  interpolation  was  used,  the  results  were  no  better.  The 
difficulty  appears  to  be  inherent  in  the  physics  of  the  situation, 
i.e.,  the  lower  the  satellite  orbit  the  more  it  deviates  from 
the  Keplerian  approximation  probably  because  of  the  non-spherical 
Earth  perturbations. 
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Variation  in  the  number  of  significant  figures  of  the  inf ut 
data  has  shown  that  5  place  accuracy  in  ephemeris  time  and  hunter 
inertial  coordinates  is  sufficient  for  most  preliminary  orbit 
determinations.  Angular  errors  of  the  order  of  +0.01°  (36  arc 

seconds)*  seem  to  be  acceptable,  but  hunter  satellite  wobble 
might  destroy  this  accuracy  unless  some  way  is  found  to  make  the 
necessary  corrections.  One  way  would  be  comparison  with  the 
local  star  field.  Errors  of  the  order  of  1°  in  either  angle 
result  in  complete  failure  of  the  method. 

The  following  table  summarizes  the  results  in  each  of  the 
orbital  configurations.  About.  6  sets  of  observations  apply  in 
each  of  the  cases  averaged. 

Given  Orbital  Element  for  Hunted  Satellite _ 


a 

e 

i 

n 

03 

M 

(a) 

10,541 

10  "  3 

10° 

0° 

90° 

135° 

(b) 

42,164 

0 

0° 

•f 

0° 

(c) 

126,500 

io~3 

10° 

0° 

90° 

0° 

Average 

Errors 

in 

Calculations 

.'a 

Ae 

Ai 

Aft 

A  w 

AM 

(a) 

3220 

0.13 

0.33° 

181° 

95° 

127 

(b) 

3300 

0.  10 

0.39° 

-20 

(c) 

6870 

0.086 

0.044 

0 

0.91° 

19.1° 

-0. 

*  This  accuracy  could  come  from  smoothing  several  observations. 
Indeterminate 
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As  suggested  above,  the  errors  in  the  computed  elements 
decrease  with  altitude,  configuration  (a)  to  (c) ,  as  the 
assumption  of  Keplerian  motion  becomes  more  accurate  (at  least 
out  to  3  times  synchronous  altitude) . 

There  is  a  so  called  neutral  point  on  a  line  connecting 
the  Earth-Moon  centers  where  the  two  gravitational  forces  are 
equal  and  opposite.  This  distance  is  approximately  345,000  km 
from  the  Earth;  as  one  approaches  this  distance, the  lunar 
perturbations  become  dominant  and  the  satellite  motion  again 
tends  to  become  non-Keplerian . 

In  configuration  (a)  the  errors  in  the  last  three  elements 
are  very  large.  The  error  in  i  is  not  bad  but  is  still  large 
enough  to  indicate  difficulty.  One  would  have  great  difficulty 
in  locating  such  a  satellite  at  some  later  time  because  of  the 
large  mean  anomaly  error.  We  find  configuration  (b)  to  yield 
more  accurate  results  even  though  some  of  the  numbers  don't 
look  much  better  than  in  (a) .  Actually  this  is  a  good  test 
because  e  and  i  are  both  zero;  therefore  ft  and  uj  are  not  defined. 
The  computed  e  and  i,  however,  are  not  zero  (hence  the  computed 
ft  and  w  will  not  be  zero) .  But  with  the  computed  inclination 
less  than  1/2°  from  the  celestial  equator  one  should  have  little 
difficulty  in  locating  the  object  again  with  a  6°  f ield-of-view 
instrument  (a  reasonably  conservative  value)  at  a  later  time 
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even  with  a  mean  anomaly  error  of  20°.  In  configuration  (c) 
all  the  errors  are  very  much  smaller  than  in  either  of  the  other 
configurations  and  very  good  preliminary  orbit  determination  is 
realized  in  this  case. 

One  example  which  has  not  been  successfully  evaluated  is 
the  highly  eccentric  orbit.  It  is  not  likely  that  one  will 
encounter  surreptitious  satellites  in  such  orbits  because  of 
the  high  risk  of  detection  at  perigee.  Nevertheless,  it  is  a 
possibility.  We  attempted  to  test  the  method  on  such  an  orbit 
but  were  unable  to  compute  the  necessary  ephemerides  because  of 
the  excessive  computer  time  needed  for  the  integrations.  We 
did  not  investigate  other  possibly  interesting  orbit  configur¬ 
ations  because  of  the  limited  time  alloted  to  this  study. 

We  conclude  that  "reasonably  accurate"  determination  of 
the  orbit  of  one  satellite  from  anotfer  by  Gauss'  method  is 
feasible  and  practicable  in  most  cases  of  interest  with  5  place 
accuracy  in  the  inertial  coordinates  and  ephernis  time  and 
+0.01°  accuracy  in  the  look  angles.  Gauss'  method  is  progressive 
more  successful  as  the  hunted  satellite  altitude  is  increased 
from  1/4  to  3  times  synchronous.  The  problem,  however,  is  a 
very  complex  one,  and  it  is  difficult  to  extrapolate  these 
results  to  all  cases  of  interest. 
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APPENDIX 


Cio_u  •  s'  Constant 

Kepler's  third  law  for  a  planet  of  mass  m  revolving  about 
the  sun  of  mass  11  may  be  stated  as  follows: 

s  0  2  1 

k" (m  +  M)T  =  4 1  a 

whore  T  is  the  period  and  a  the  semi-major  axis.  The  value  of 
k  depends  upon  the  units  of  time,  distance  and  mass.  Following 
Gauss,  the  units  are  taken  to  be  the  solar  mass,  the  mean  solar 
day  and  the  Earth's  mean  distance  from  -the  Sun  (the  astronomical 
unit)  in  heliocentric  problems.  Kepler's  third  law  becomes 

2  2  2  1 
K  (1  +  m)T  =  4 ii  a 

The  quantity  k  is  called  the  Gaussian  constant  of  gravitation 
and  in  the  units  defined  above  we  have  k  =  0.01720209895. 

For  satellite  motion  about  the  Earth  we  can  choose  any  con¬ 
venient  units:  the  ephemeris  minute,  the  mass  and  radius  of 
the  Earth  or  the  ephemeris  day,  the  mass  of  the  Earth  and  km, 
etc.  To  determine  the  value  of  k  in  this  latter  case  we  use 

Kepler's  third  law  and  take  m  to  be  the  mass  of  the  Moon  in 

_1_ 

units  of  Earth  masses ,  (81 . 3*1)  ,  T  to  be  27.321601  days  and  a 
to  be  377  ,028  km,  the  semi -major  axis  of  tlie  lunar  orbit.  Then 
the  quantity  k  lias  a  value  5.2915  E+7.  If  we  use  the  astro¬ 
nomical.  unit  instead  of  km,  then  k  has  a  value  2.89125  E-5. 
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